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We study the rheology of confined suspensions of neutrally buoyant rigid monodisperse spheres 
in plane-Couette flow using Direct Numerical Simulations. We find that if the width of the channel 
is a (small) integer multiple of the sphere’s diameter, the spheres self-organize into two-dimensional 
layers that slide on each other and the suspension’s effective viscosity is significantly reduced. Each 
two-dimensional layer is found to be structurally liquid-like but their dynamics is frozen in time. 


Suspensions of solid objects in simple Newtonian sol¬ 
vents (e.g., water) can show a kaleidoscope of rheological 
behaviors depending on the shape, size, volume fraction 
(0) of the additives, and the shear-rate (7) imposed on 
the flow; see, e.g., Refs. [1, 2] for reviews. Suspensions 
can be of various types, e.g., suspensions of small parti¬ 
cles (smaller than the viscous scale of the solvent), where 
the solvent plays the role of a thermal bath, are called 
Brownian suspensions (e.g. colloids) [1]. At small <f> 
and under small 7, the effective viscosity of such sus¬ 
pensions increases with <f>: p e ff ~ (5/2)0 [3], as de¬ 

rived by Einstein [4] (See also Ref. [5] for a d dimen¬ 
sional generalization). In the other category are suspen¬ 
sions of large particles (e.g. emulsion, granular fluids 
etc.), where the thermal fluctuations are often negligi¬ 
ble. Such suspensions are called non-Brownian suspen¬ 
sions. For moderate values of 0 and 7 non-Brownian 
suspensions show continuous shear-thickening , (i.e., /z e ff 
increases with 7) [6, 7] which can be understood [7] in 
terms of excluded volume effects. Such a rheological re¬ 
sponse has been observed in many natural and industrial 
flows, including flows of mud, lava, cement, and corn¬ 
starch solutions. Dense (large 0, close to random close 
packing) non-Brownian suspensions may show discontin¬ 
uous shear-thickening [8, 9] - a jump in effective viscosity 
as a function of 7. 

Recent experiments have uncovered intriguing rheo¬ 
logical behavior of very dense (large 0) suspensions un¬ 
der confinement [6, 8, 10]. Common wisdom dictates 
that, under confinement, the inertial effects are gener¬ 
ally unimportant at small 7. But a series of recent stud¬ 
ies [11-14] have demonstrated that the effect of fluid in¬ 
ertia, although small, can give rise to variety of effects 
even in microfluidic flows. Drawing inspirations from 
these two sets of works, in this paper, we study the ef¬ 
fects of confinement on non-Brownian suspensions with 
moderate values of 0 and 7. In particular, we choose the 
range of 0 and 7 such that in bulk the suspensions show 
continuous shear-thickening. 

As a concrete example, we consider direct numerical 
simulations (DNS) of three-dimensional plane-Couette 
flow - with the x, y and z coordinates along the stream- 


wise, span-wise and wall-normal directions respectively 
[see Fig. (la)] - embedded with neutrally-buoyant rigid 
spheres of radius a. The fluid phase is described by solv¬ 
ing the incompressible Navier-Stokes equations in three 
dimensions. The fluid is sheared by imposing constant 
stream-wise velocity of opposite signs, Uo = 7L2/2 at 
2 = ±L Z / 2. Periodic boundary conditions are imposed 
on the other two directions (x and y with lengths L x and 
L y respectively). The motion of the rigid spheres and 
their interaction with the flow is fully resolved by using 
the Immersed Boundary Method [15, 16]. A description 
of the equations and the details of the algorithm is pro¬ 
vided as supplemental material. 

We study the effects of confinement by changing the 
dimensionless ratio £ = L z /2a, where L z is the channel 
width in the 2 direction. In practice, we change L z but 
hold the particle radius a fixed. The effective viscosity, 
is thus function of the dimensionless numbers, 0, £ 
and the particle Reynolds number, Re = /rya 2 //if, where 
p and /if are the density and dynamic viscosity of the sol¬ 
vent. The most striking result of our simulations is that 
at or near integer values of £, /i e ff decreases significantly 
compared to its bulk value; see Fig. (lb). This drop 
can be so large that the net rate of mass transport for a 
thinner channel (£ = 2) is more than in a wider channel 
(£ = 2.5), see Fig. (lc). We further demonstrate that, at 
small integer values of £, the rigid spheres self-organize 
into an integer number of horizontal layers, which slide 
on one another, see Fig. (2). This, in turn, decreases 
the transport of horizontal momentum across the layers 
generating the drop in effective viscosity. We also show 
that at integer values of £, for which layers appear: (a) 
the movement of spheres across layers is not a diffusive 
process, Fig. (3), (b) the typical residence time of the 
spheres in a layer is much longer than the time scale set 
by the shear, 7 -1 ; in fact, a large number of spheres 
never leave the layer within the runtime of the simula¬ 
tions Fig. (4a); (c) the layers are structurally liquid-like 
but dynamically very slow compared to the timescale set 
by 7 -1 , see Fig. (4). 

Initially, the spheres are placed at random locations, 
with no overlap between each other, and with velocities 
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FIG. 1. (a) A sketch of the domain, (b) Surface plot of the effective viscosity f as a function of £ and Re for £> = 0.3. 

The inset shows /x e ff//^f versus £ for £> = 0.3 and for three different values of Re = 1 (red circles), 5 (blue squares) and 10(black 
triangles), (c) The flow rate of matter (both solvent and additive ) F (blue squares, left vertical axis), defined in (1), and the 
flow-rate-per-unit-cross-section, /, (red circles, right vertical axis) as a function of £, for Re = 5 and £> = 0.3. 


that are equal to the local fluid velocity, which is taken 
to be the laminar Couette profile, Fig. (la). We calculate 
the effective viscosity, /i e fr, as the ratio between the tan¬ 
gential stress at the walls and the shear 7 averaged over 
the walls and over time. See supplemental material for 
further details of measurement and estimate of errors. 

The effective viscosity is shown in Fig. (l)b as 
a function of Re and £, for £> = 0.3. For large chan¬ 
nel widths, e.g., £ = 5, and large volume-fractions, e.g, 
£> > 0 . 2 , we obtain shear-thickening [see Supplemental 
Material Fig. ( 6 )] as was found earlier in Ref. [7]. 

Here we address how confinement affects the rheol¬ 
ogy. Experiments [17], in agreement with earlier analyt¬ 
ical calculations [18], have found that at small volume- 
fraction ( (j) < 0.15), /i e ff increases as £ decreases. Our 
results [Supplemental Material, Fig. ( 6 )] also supports 
this conclusion. Furthermore, our simulations can ac¬ 
cess hitherto unstudied higher values of volume-fraction 
(£> > 0 . 2 ) for which this trend seems to reverse, i.e., /i e ff 
decrease with confinement, see e.g., Supplemental Mate¬ 
rial, Fig. ( 6 ). On deeper scrutiny, a more striking picture 
emerges. As we show in Fig. (lb), for £> = 0.3, at or near 
small integer values of £ (£ = 2,3, and 4), the effective 
viscosity drops significantly. In our simulations, the min¬ 
imum value of viscosity is found at £ = 2 , which is as low 
as 50% of its bulk value (£ > 6 ). 

To appreciate how dramatic this effect is, we measure 
the (dimensionless) flux of matter (i.e. both the fluid and 
the spheres) through the channel, defined as [19] - 

F = a 2 p{ J U x p(z)dydz); ( 1 ) 

where U = £Z7 P + (1 — £)u, with £ = 1 at a grid point 
inside a rigid sphere but £ = 0 otherwise [20], Here, U p 
is the velocity of the sphere, u is the velocity of the fluid, 
p(z) = 1 for z > 0 and — 1 for z < 0 , and (•) denotes aver¬ 
aging over time. In Fig. (lc) we report F as a function of 
£. As the flux F is not normalized by the cross-sectional 


area of the channel it is expected to increase linearly as 
a function of £. This expectation indeed holds for £ > 4. 
But, below that, for integer values of £ the effective vis¬ 
cosity can decrease so much that F is not even a mono¬ 
tonic function of £, in particular, F(£ = 2) > F(£ = 2.5); 
the flux through a wider channel is actually smaller. The 
flux-per-unit-area, f = F(a 2 /L y L z ), shown in Fig. (lc) 
is expectedly a constant at large £. For the small integer 
values of £, it is significantly higher than its bulk value. 

To investigate the mechanism behind the rheology, we 
examine snapshots of the spheres, see Fig. (2a) for £ = 2 
(top) and £ = 2.5 (bottom). The spheres are color-coded 
by their initial wall-normal locations (red corresponds to 
an initial position near the top boundary and blue to the 
bottom boundary). It can be clearly observed, that for 
£ = 2, the particles form a bi-layered structure. This lay¬ 
ering is also confirmed by the wall-normal profiles of the 
average particle number density displayed in Fig. (2b) for 
£ = 2, 2.5 and 3. In the first and the last case, one can 
observe equally-spaced two and three prominent peaks 
respectively. Note that a weaker layering is observed 
for £ = 2.5. The drop in effective viscosity thus corre¬ 
sponds to the formation of layers that slide on each other, 
with little transport of momentum across the layers. For 
£ = 2 , where layering is most prominent, the particles 
form disordered liquid-like structures, within each layer, 
as seen by the radial distribution function of the position 
of the spheres, see Supplemental Material, Fig. (7). 

In order to understand the dynamics of particles in the 
wall-normal direction, we show in Fig. (3a), the proba¬ 
bility distribution function (PDF) of the displacement of 
the spheres, d z (t) = z c (t)—z c (0 ), at different times, where 
z c (t) is the the z coordinate of the center of a sphere. Ob¬ 
viously, as t -A 0 the PDF [P(d z )\ must approach a Dirac 
delta function. Remarkably, for £ = 2 the PDF has expo¬ 
nential tails, i.e., it is non-diffusive, with some particles 
undergoing larger displacements as shown in Fig. (3a). 
At later times, for £ = 2 , the PDF of d z develops a 
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FIG. 2. (a) Three dimensional view of positions of the spheres 
for Re = 5, £> = 0.3 and £ = 2 (top) and £ = 2.5 (bottom) at 
later times after the initial transients. The particles are color 
coded by their initial wall-normal position, (red: close to 
the top boundary, blue: close to the bottom boundary), (b) 
Average local density of spheres p p versus the wall-normal 
coordinate z, for £ = 2 (red triangles), 2.5 (blue squares), and 
3 (black dots). 


peak at d z /L z = 1, indicating the hopping between two 
layers. A similar behavior is observed for £ = 3 too. 
Contrast this result with the PDF of the displacement 
along the span-wis e(y) direction, inset of Fig. (3a), which 
clearly shows Gaussian behavior at all times, implying 
diffusive dynamics. The second moment of these PDFs 
provides the mean squared displacement of the particles, 
S^it) = {[z c (t) — £ c ( 0 )] 2 )p> the time evolution of which is 
shown in Fig. (3b). At late times, in general, a power-law 
dependence on time is found, S^if) ~ t 13 , where /3 = 1 
would imply a simple diffusive behavior. For £ = 2, 3 we 
estimate /3 « 0.61 and 0.88 respectively. But for £ = 2.5 
, ~ 1 is obtained. The diffusive behaviour for £ = 2.5 

can be further confirmed by plotting the PDF of d z . At 
intermediate times (t = 2.57 -1 ), for £ = 2.5, the PDF 
develops Gaussian tails, indicative of a diffusive process 
and at late times it approaches a constant (see Supple¬ 
mental Material, Fig. (8)). 

To visualize the dynamics, we provide movies 
of the particles’ trajectories (available at https:// 
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FIG. 3. (a) The PDF, P(d z ), of the displacement (d z ) of 

the center of the spheres along the z direction for £ = 2, at 
early times, t = 2.5'y~ 1 (red, filled triangles) and late times 
t = 386 1 (red, open triangles). For comparison, the inset 
shows the PDF, P(d y ), of the displacement (d y ) along the y 
direction at early times, t = 2.5y _1 (black, filled diamonds) 
and late times t = 5y -1 (black, open diamonds) (b) Log-log 
plot of mean square displacement of the spheres, S'!, versus 
t r y. At late times, S| ~ where /3 = 1 implies diffusion. 
For three different values of £ = 2 (red triangles), 2.5 (blue 
squares) and 3 (black dots), we obtain /? « 0.61, 0.95 and 
0.88, respectively. 


www.youtube.com/wat ch?v=Qn4DiXZFsbw for the case 
£ = 2 and at https://www.youtube.com/watch?v= 
AmNsAsY0eC8 for the case £ = 2.5 ). These clearly demon¬ 
strate that for £ = 2 , each horizontal layer is structurally 
disordered but dynamically frozen. We quantify this phe¬ 
nomenon by three different measurements: 

(A) We calculate the residence time (r) of a sphere in 
a single layer. A sphere is considered to reside within 
a horizontal layer till the wall-normal coordinate of its 
center, z c , is within a distance of 2 a of its initial position. 
Instead of calculating the PDF by histograms of r, we 
calculate the cumulative PDF, Q(r), by the rank-order 
method [ 21 ], as the latter is free from binning errors. 
The cumulative PDF, Q(r), is displayed in Fig. (4a) as 
a function of r for £ = 2 and 2.5 for Re = 1 and 5. For 
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FIG. 4. (a) Cumulative PDF, Q, of the residence times of spheres (yr) in a layer for £> — 0.3, Re = 1 (in blue) , and Re = 5 
(in black) and £ = 2 (circles), and £ = 2.5 (squares), (b) The streamwise velocity-velocity auto-correlation function, R xx , for 
Re = 5 and <j) = 0.3 for four different values of £ = 2(red triangles), 2.5 (blue squares), and 3 (black dots), (b) Evolution of 
averaged square of relative distance between pairs of sphere, (rjj(£, ro)) p which were at distance ro at t = 0; for ro = 4a, for 
three different values of £ = 2 (red triangles), 2.5 (blue squares), and 3 (black dots). 


£ = 2, Q remains very close to unity during the whole 
duration of the simulation, i.e., very few spheres actually 
move from one layer to another. In other words, the 
layers are quite stable to perturbations in wall-normal 
directions. Conversely for £ = 2.5, Q(r) can be fitted by 
a Gaussian. 

(B) The streamwise velocity auto-correlation function 
of the spheres, R xx = (U£ (t) (0)), is shown in Fig. (4)b. 
For £ = 2, R xx ^ 1 for a very long time, implying that the 
temporal fluctuations of the stream-wise velocity are neg¬ 
ligible. This suggests that each sphere moves in a layer 
with a uniform stream-wise velocity keeping their rela¬ 
tive distances practically constant. For the cases where 
the layering is not very strong, e.g., £ = 2.5 the auto¬ 
correlation function decays rapidly. For £ = 3, layering 
reappears and R xx again shows slow decay in time. 

(C) Let us define ry(t;ro) to be the (horizontal) dis¬ 
tance between a pair of spheres at time £, which were at 
a (horizontal) distance ro and t = 0. If the layers formed 
by the spheres were truly frozen-in-time we would ob¬ 
tain r||(£;ro) = ro for all t and ro- In Fig. (4c), we show 
the time evolution of (rjj(t; ro)) p , for ro = 4a, where the 
symbol (*) p denotes averaging over all possible particle 
pairs [22]. Had the spheres moved chaotically within 
a layer, ry would have grown exponentially with time. 
Clearly for all the cases shown in Fig. (4c), ry grows at 
most linearly with time. In particular, when layering 
occurs (£ = 2,3), it takes a long time (t > 4307 -1 ) for 
ry (£; 4a) to grow by a factor of two. This quantifies again 
the dynamical stability of the layer, the relative in-plane 
distance between pairs of spheres change slowly (linearly) 
with time. 

In conclusion, using numerical simulations, we demon¬ 
strate that the effective viscosity of an extremely confined 
non-Brownian suspension can exhibit a non-monotonic 
dependence on the channel width, in particular the effec¬ 
tive viscosity sharply decreases if the channel width is an 
(small) integer multiple of particle diameter. We demon¬ 


strate that this behavior is accompanied by a change in 
micro-structure, namely the formation of particle layers 
parallel to the confining plates. The two-dimensional lay¬ 
ers formed by the particles slide on each other, the layers 
are structurally liquid-like by evolve on very show time 
scales. Similar layering under shear, have been theoreti¬ 
cally anticipated [23], but the consequences for transport 
in extreme confinement as shown by us has never been 
demonstrated before. We finally note that our results are 
in contrast with the case of sheared Brownian suspensions 
where the structure was found to be uncorrelated with 
measured viscosity [24]. 
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Supplemental Material 

Here we provide supplemental material for our paper. 
The material is divided into several subsections. 

Direct numerical simulation 

Simulations have been performed using the numerical 
code originally described in Ref. [25] that fully describes 
the coupling between the solid and fluid phases. The 
Eulerian fluid phase is evolved according to the incom¬ 
pressible Navier-Stokes equations, 

V • u = 0 (A.2) 

d t u + u • V'u = — - Vp + v\/ 2 u + / (A.3) 

P 

where u, p and v = p/p are the fluid velocity, density 
and kinematic viscosity respectively {p is the dynamic 
viscosity), while p and f are the pressure and a generic 
force field (used to model the presence of particles). The 
particles centroid linear and angular velocities, U p and 
f2 p are instead governed by the Newton-Euler Lagrangian 
equations, 

d/7 p 7 

Pp v p~JT = P f T-ndS (A.4) 

at Jdv p 

^ dOP r 

2 p—— = p (p r x r • ndS (A.5) 

dt Jdv p 

where V p = 47ra 3 /3 and X p = 2p p E p a 2 /5 are the parti¬ 
cle volume and moment of inertia; r = —pi + 2 pE is 
the fluid stress, with E = (Vn + Vu T ) /2 the defor¬ 
mation tensor; r is the distance vector from the center 
of the sphere while n is the unity vector normal to the 
particle surface dV p . Dirichlet boundary conditions for 
the fluid phase are enforced on the particle surfaces as 
u\qv p = U p + x r. 

An immersed boundary method is used to couple the fluid 
and solid phases. The boundary condition at the moving 
particle surface is modeled by adding a force field on the 
right-hand side of the Navier-Stokes equations. The fluid 
phase is therefore evolved in the whole computational do¬ 
main using a second order finite difference scheme on a 
staggered mesh while the time integration is performed 
by a third order Runge-Kutta scheme combined with a 
pressure-correction method at each sub-step. This in¬ 
tegration scheme is also used for the Lagrangian evolu¬ 
tion of eqs. (A.4) and (A.5). Rach particle surface is de¬ 
scribed by uniformly distributed TVl Lagrangian points. 
The force exchanged by the fluid on the particles is im¬ 
posed on each l — th Lagrangian point. This force is 
related to the Eulerian force field f by the expression 
fix) = Yld =i — Xi) A Vi (where AV\ is the volume 
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of the cell containing the l — th Lagrangian point while S 3 
is the three dimensional Dirac delta distribution). The 
force field is calculated through an iterative algorithm 
that ensures a second order global accuracy in space. In 
order to maintain accuracy, eqs. (A.4) and (A.5) are 
rearranged in terms of the IBM force field, 

PpVp = -p F ' 1 + Pj t f u dv (A - 6) 

1=1 ^ y p 

dftP A „ ATr d f 

^ =-' , 5> x F ' AVt+ fj t L r x udV 


where r\ is the distance from the center of a particle. The 
second terms on the right-hand sides are corrections to 
account for the inertia of the fictitious fluid contained in¬ 
side the particles. Lubrication models based on Brenner’s 
asymptotic solution gare used to correctly reproduce the 
interaction between two particles when the gap distance 
between the two is smaller than twice the mesh size. A 
soft-sphere collision model is used to account for colli¬ 
sions between particles and between particles and walls. 
An almost elastic rebound is ensured with a restitution 
coefficient set at 0.97. These lubrication and collision 
forces are added to the right-hand side of eq. (A.6). 

The code has been validated against several classic test 
cases [25] and has been used earlier to study shear¬ 
thickening in inertial suspensions in Ref. [7]. 


Measurement of effective viscosity 

We calculate the effective viscosity as the ratio between 
the tangential stress at the walls and 7 . The tangential 
stress, and consequently the effective viscosity, is differ¬ 
ent at different points of the wall and also at each point 
changes as a function of time. At each time, we calcu¬ 
late the average of this effective viscosity over the walls 
and obtain a time-series fjL e s(t). This time-series is dis¬ 
played in Fig. (5) for one particular case. Clearly, after 
a short while the effective viscosity reaches a stationary 
value with fluctuations about it. Typically, we find that 
our simulations reach statistically stationary state when 
time t > 2007 “ 1 The average of yU e ff(^) over this statis¬ 
tically stationary state is the effective viscosity /i e ff and 
the standard deviation of the fluctuations is used as an 
estimate of the error in the measurement of fi e ^ as re¬ 
ported in Tablel A careful look at the table will convince 
the reader that the relative strength of the fluctuations 
decreases at commensuration. 


£ (Re = 1 ) Meff/Mf (Re — 5) (R e — 10) 


1.5000 

4.665 ±0.113 

5.123 ±0.079 

5.161 ±0.107 

1.7500 

5.758 ±0.116 

5.925 ±0.114 

6.100 ±0.115 

1.8750 

5.334 ±0.186 

5.603 ±0.190 

5.895 ±0.173 

1.9375 

4.641 ±0.190 

4.906 ±0.234 

5.213 ±0.209 

2.0000 

3.650 ±0.120 

3.709 ± 0.058 

3.831 ±0.066 

2.0625 

3.888 ±0.152 

3.916 ±0.143 

3.596 ±0.075 

2.1250 

4.237 ±0.183 

3.928 ±0.139 

3.574 ±0.100 

2.2500 

4.517 ±0.141 

4.308 ± 0.244 

3.581 ±0.087 

2.5000 

5.129 ±0.165 

5.396 ±0.189 

5.836 ±0.181 

2.7500 

4.403 ±0.162 

5.083 ±0.163 

5.766 ±0.158 

2.8750 

4.381 ±0.147 

4.609 ±0.138 

5.051 ±0.148 

2.9375 

4.303 ±0.147 

4.261 ±0.120 

4.223 ±0.102 

3.0000 

4.506 ±0.155 

4.446 ±0.130 

4.148 ±0.180 

3.0625 

4.462 ±0.177 

4.763 ±0.148 

4.326 ±0.223 

3.1250 

4.718 ±0.185 

4.942 ±0.171 

4.670 ± 0.237 

3.2500 

4.721 ±0.184 

5.286 ±0.160 

5.449 ±0.172 

3.5000 

4.608 ±0.165 

5.535 ±0.166 

6.238 ±0.189 

3.7500 

4.610 ±0.168 

5.326 ±0.142 

5.923 ±0.166 

4.0000 

4.629 ±0.156 

5.323 ±0.154 

5.358 ±0.176 

5.0000 

5.150 ±0.208 

5.626 ±0.152 

6.371 ±0.193 

6.0000 

4.752 ±0.180 

5.816 ±0.143 

6.807 ±0.173 


TABLE I. Values of obtained at different £ and Re 

with respective error. 


(a) 



(b) 



FIG. 5. The effective viscosity fi e ^{t) (normalized by the vis¬ 
cosity of the solvent) as a function of time for one represen¬ 
tative run. 
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FIG. 6. The effective viscosity fies/fi f versus 0 for two differ¬ 
ent values of Re = 1 (dashed lines) and 5 (continuous lines) 
for three different values of £ = 1 (black filled circles), 2 (blue 
squares) and 5 (red triangles). 


Effective viscosity as a function of volume fraction 

Here we show, see Fig. (6), the variation of the effective 
viscosity, /i e ff as a function of the volume fraction, <fi. 

In-layer radial distribution function 

The radial distribution function of the position of the 
centers of the spheres in the x — y plane is shown in 
Fig. (7). This is defined as 

_ 1 dN r . 

R \| = ——---— (A.8) 

11 27rrAzn 0 dr 

where N r is the number of particles is the number of 
particles in a cylinder of radius r and no = N p (N p — 
1)/(2V) is the density of particles pairs in a volume V 
with N p the total number of particles. 
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■£=2.0, Re=5 


FIG. 7. The radial distribution function of the position of the 
centers of the spheres in the x — y plane. 



FIG. 8. The PDF (P(d z )) of the displacement (d z ) of the 
center of the spheres along the z direction, for £ = 2.5, at 
early times, t — 2.5y _1 (red, filled triangles) and late times 
t = 3867 “ 1 (red, open triangles). For comparison, the inset 
shows the PDF (P(d y )) of the displacement (d y ) along the y 
direction at early times, t = 2.5y _1 (black, filled diamonds) 
and late times t = 57 “ 1 (black, open diamonds) 
















